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Dynamo coefficients from local simulations of the turbulent ISM 
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Observations in polarized emission reveal the existence of large-scale coherent magnetic fields in a wide range of spiral 
galaxies. Radio-polarization data show that these fields are strongly inclined towards the radial direction, with pitch angles 
up to 35° and thus cannot be explained by differential rotation alone. Global dynamo models describe the generation of 
the radial magnetic field from the underlying turbulence via the so called a-effect. However, these global models still rely 
on crude assumptions about the small-scale turbulence. To overcome these restrictions we perform fully dynamical MHD 
simulations of interstellar turbulence driven by supernova explosions. From our simulations we extract profiles of the 
contributing diagonal elements of the dynamo a-tensor as functions of galactic height. We also measure the coefficients 
describing vertical pumping and find that the ratio 7 between these two effects has been overestimated in earlier analytical 
work, where dynamo action seemed impossible. In contradiction to these models based on isolated remnants we always 
find the pumping to be directed inward. In addition we observe that 7 depends on whether clustering in terms of super- 
bubbles is taken into account. Finally, we apply a test field method to derive a quantitative measure of the turbulent 
magnetic diffusivity which we determine to be ~ 2 kpckms -1 . 
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1 Introduction 

The last decade with its advances in observational 
technology has brought increasingly detailed maps of 
galactic magnetic fields. Radio-polarization data indicate 
large scale regular fields within the interstellar medium. 
For typical spiral ga l axies the ISM is highly turbulent 
<|Mac Low & KlessenI |2004 . The main drivers of the 
turbulence are thought to be wi nds of massive stars and ex- 
plosions of supernovae (SNe) ( de Avillez & Breitschwerdtl 



2005; Uoung & Mac Low! 120061) . as well as 



active region s) the magneto-rotational 
ity (MRI) ( Dziourkevitc h. Elstner & Ri idiger 



(in less 
instabil- 
2004 



Piontek & Ostriker 20051, 120071) . Furthermore. lHanasz et al, 
d2004 l2006l) ~ have shown that a Parker-type instability 
can be driven by cosmic rays, an idea first advocated by 
Parker I d 1992b . It, however, remains a matter of discussion, 
whether c osmic rays play a sig nificant role in the interstellar 



wnetner cosmic rays play a sigi 
dynamics (ISnodin et al . 2006). 



The generation of a mean magnetic field from turbu- 
lent fluctuations can be explained via the so called a-effect 
which describes the correlations of the small-scale turbu- 
lent velocity and magnetic field giving rise to a mean elec- 
tromotive force (EMF). In the case of the ISM there are 
three characteristic asymmetries leading to non-vanishing 
correlations: (i) the axis of rotation, (ii) the galactic shear 
gradient, and (iii) the (vertical) gradient in den sity and tur- 
bulence intensity (Riidiger & Kichatin ovll993 ). The mutual 
strength of the different terms puts constraints on the oper- 



ability of a dynamo process. The ratio 7 of the diamagnetic 
pumping over the a-effect has an influence on the efficiency 
of the dynamo. The strength of the a-effect compared to the 
shear puts limits on the pitch angle. 

Until recently, a direct numerical simulation of the tur- 
bulent ISM has been infeasible. The refore, early theoretical 
model s like the SOCA approach by Riidig er" fe Kichatinovl 
( 11993b predicted the outcome of ISM-turbulence based 
on simplifying assumptions. In an alternative approach 
iFerrierd {1992) derived the dynamo effect for isolated su- 
pernova remnants (SNRs) and super-bubbles (SBs). How- 
ever, these early no-interaction models arrived at pro- 
hibitively high values for 7. To test these findings first nu- 
merical si mulations for single SNR s have been performed 
in 2D bv iKaisig. Riidiger & Yorkd (1 19931) and in 3D by 
Ziegler. Yorke & Kaisigl (119960 . Still the key issue of a 



dominating turbulent pumping remained. The situation was 
somewhat imp roved by taking into account stratification 
dFerriere| [l998. hereafter FER98) yielding a value 7 « 6. 
The parameter ra nge for 7 permitting dynamo soluti ons has 
been explored bv lSchultz. Elstner & Riidigerl dl994 . 



The major limitation to the no-interaction models de- 
scribed above results from the fact that (even for SBs) the 
explosions are considered as isolated events taking place on 
a uniform background. But this is not the c ase for the ISM 
which due to thermal instability ( Field! 19651) is a highly het- 
erogeneous medium. To overcome these limitations we per- 
form direct numerical simulations of the local ISM and use 
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a test field method ( Sch rinner et al.l 120 05. 2007) to obtain 
the dynamo parameters. 
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2 Physical model and parameters 

We simulate the dynamic evolution of the stratified, turbu- 
lent ISM utilizing a 3D MHD model including the physical 
effects described below. The computational domain covers 
a box of 0. 8x0. 8x4. kpc 3 vertically centered around the 
midplane and representing a local patch of the galactic disk. 
The purely vertical gravi tational potential is adopted from 
Kuijken & Gilmord (Il989h . 

The equations of resistive MHD are solved in the local 
shearingbox approach, i.e., we apply a co-rotating Carte- 
sian coordinate system with x, y, and z being the unit vec- 
tors along the radial, azimuthal, and vertical direction. The 
background shear of the flow is characterized by the pa- 
rameter q = dlnfl/dhiR, where the case q = — 1 cor- 
responds to a flat rotation profile. Additional source terms 
Tsn — P 2 A(T) + pT(z) in the total-energy equation repre- 
sent the thermal energy input due to supernovae and opti- 
cally thin radiative cooling/heating. 

2.1 Supernova driving 

We model supernova explosions as local injections of ther- 
mal energy resulting in turbulence at mildly supersonic 
Mach numbers. To assure convergence of the emerging su- 
pernova remnants numerical solutions with increasing spa- 
tial resolutio n As have been tested against the ana lytical de- 
scription bv lCioffi. McKee & Bertschingerl (1198a) . The re- 
sults agree wel l and c orrespond to the findings reported in 
Mac Low etal] d2005l) . 

SN-events are exponentially distributed in the verti- 
cal direction with scale heights of 325 pc for type-I and 
90 pc for type-II SNe. For the models presented in this pa- 
per we use 3/4 the galactic frequencies which are o\ = 
4Myr~ 1 kpc~ 2 respectively on = 30 Myr -1 kpc -2 . The 
associated exp losion energies are 10 51 and 1.14 x 10 51 erg 
dFerrierell200ll) . 

Within our model we make an important distinction be- 
tween type-I and type-II SNe. The latter are spatially clus- 
tered by the (artificial) constraint that the density at the 
explosion site must be above average (with respect to a 
horizontal slab) wh ile the former are spatially uncorrelated 



dKorpi et al. 19991) . We use this simple prescription as a 



prox y for a more self-consistent treat ment of the cluster- 
ing ( de Avillez & Breitschwerdtl 2005 ) and find a fraction 
of clu stered events comparable to observations dFerriere 
2001). The reference simulation without SBs reveals that 
the general morphology is affected quite strongly by the 
clustering indicating the importance of this effect. 

2.2 Radiative cooling and diffuse heating 

We treat the interstellar medium as an optically thin 
plasma and prescribe the coupling to the radiation field 
via a piecewise power law of the form: A(T) = A^T^, 
for Ti < T < Ti + \. The parameters used are es- 
sentially a combination of the cooling curves given by 
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Table 1 Overview of conducted models. Clustering ('cl.') 
applies to type-II SNe only. SN-rates are 3/4 the galactic 
values. 



Sanchez-Salce do. Vazquez-Semadeni & Gazol (2002) and 
Slvz et alj ( 12005b . where the coefficients have been slightly 
modified to make th e resulting curve co ntinuous. In contrast 
to previous work by Korpi et al. I dl999h we include the ther- 
mally unstable regime below 6102 K leading to the forma- 
tion of a cold ISM phase. To numerically resolve the cooling 
instability we apply thermal conduction such that the Field 



length Af is covered by 4 grid cells (cf. iKovama & Inutsuka 



2004). The various effects that contribute to the diffuse heat- 
ing of the interstellar plasma are subsumed in a prescribed 
heating r ate T(z) = Tn e~ z / h , wi th h = 300 pc as dis- 



cussed in lJoung & Mac Low! (12006). 



2.3 The initial model 

Previous stratified models 



de Avillez & Breitschwerdt 



(Korpi etal. 1999 



2004 Uoung & Mac Low 



2006; Pio ntek & Os triker 2007) all start from an isothermal 
initial state at a prescribed temperature. The main drawback 
of this is that the isothermal stratification is not in radiative 
equilibrium and the disk will instantaneously collapse 
until the dynamic pressure from SN- or MRI-turbulence 
will balance this process. To avoid this undesired behavior 
we propose a more sophisticated initial model where the 
vertical profiles of density and pressure are numerically 
integrated to be in combined hydrostatic and radiative 
equilibrium. The computed profiles are considerably flatter 
than the isothermal ones while the temperature varies by a 
factor of about five. 

For our fiducial model we choose an initial midplane 
density po corresponding to one particle per cm 3 which re- 
sults in an equilibrium pressure of po/^s = 6000 K cm -3 . 
We use a mean molecular weight of p, — 0.6 assuming a 
fully ionized plasma of cosmic abundance. For the current 
study the rotation rate is fixed at f2 = 1 00 km s ~~ 1 kpc ~ 1 . 
As in iBrandenburg. Korpi & Mee (2007) we scale the dy- 
namic viscosity coefficient with the density, i.e., we use a 
constant kinematic viscosity of v = 5 x 10 24 cm 2 s _1 . To 
obtain a constant Prandtl number of Pr = 4.2 we apply the 
same scaling to the thermal conduction coefficient k. 

The magnetic Prandtl number Pm = vjr\ is thought to 
be very high for the ISM. With the limited dynamic range of 
our simulations we are, however, restricted to values close 
to unity. For practical purposes we choose Pm = 2.5 equiv- 
alent 77 = 2xl0 24 cm 2 s _1 which is still two orders of mag- 
nitude smaller than the expected turbulent diffusivities. 
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Fig. 1 Vertical slices of the top half of our box (upper panels), and horizontal slices through the midplane (lower panels) 
after t — 112 Myr. Quantities shown are: (a) number density [cm -3 ], (b) column density [cm -2 ], (c) temperature [K], 
(d) velocity dispersion [kms -1 ], and (e) magnetic field strength [/xG]. The logarithmic grey scales extend over ranges 
[-4.82, 1.10], [17.34, 21.55], [2.06, 7.20], [-0.15, 2.77], and [-4.16, 0.34] for the corresponding variables. 



The vertical boundary conditions are of the outflow 
type with the magnetic field extrapolated according to the 
solenoidal constraint. Due to the (sheared-) periodic bound- 
ary conditions in the x- and y-direction the vertical mag- 
netic field is ideally conserved for any horizontal slab. This 
means that there are two distinct classes of models includ- 
ing or excluding vertical magnetic flux. For simplicity we 
focus on the latter case and start from an initially toroidal 
plus radial field with a pitch angle of —6°. The field is 
scaled with \/p/po to yield a constant plasma parameter 
/3 P = 2p/B 2 w 2000 throughout the disk. The resulting 
Alfven velocity for the initial model is 0.3 — 0.6 km s -1 
whereas the sound speed ranges from 10 — 20 km s -1 . The 
particular choice of parameters for the conducted models is 
summarized in Table [TJ 



3 Numerical methods 

For our computations we make use of the newly dev eloped 
version 3 of the NIRVANA code JZieglerl2004l2005h which 
is a general purpose MHD fluid tool employing the tech- 
nique of adaptive mesh refinement (AMR). For the simula- 
tions presented here this feature is, however, not being em- 
ployed. We extended the standard shearingbox model to be 
compatible with the conservative numerical scheme and ap- 



ply flux-matching at the sheared interfaces to im prove the 
conservation properties ( Gressel & Zieglenl2007h . 

3.1 Dynamo test fields 

To find a closure for the mean field equations one 
strives to parameterize £ = u' xB' with respect to av- 
eraged quantities. Here we adopt the standard descrip- 
tion where £ depends on the mean field and its gradi- 
ents. When using spatia l avera ges along horizontal slabs 
Brandenburg & Sokoloffl (120021) showed that the resistivity 
tensor can be reduced and one yields: 

£i = aijBj - fjijSjkidk&i , i,j e{R,<t>},k = z. (1) 

While the diagonal elements of a. quantify the dynamo 
process, its anti-symmetric off-diagonal elements consti- 
tute the vertical pumping effect described by the parame- 
ter 7 Z = 0.5 {cLfyR — OLR^). Similarly the diagonal elements 
of f) are interpreted as turbulent resistivity r\ t , while its off- 
diago nal components can lead to fl x J-type dynamo ef- 
fects jRtidiger & Hollerbachll2004 . In total we have 4+4 
unknowns that we wish to determine from equation (Q]). 

To avoid complications with the inversion of this tenso- 
rial equation, we apply the t est field approach proposed by 
Schrinner et al. I d2005l 2007). The method has also rec ently 
been adopted to the shearingbox case by iBrandenburg 



www.an-journal.org 



© 2008 WILEY-VCH Verkg GmbH & Co. KGaA, Weinheim 



4 



O. Gressel, U. Ziegler, D. Elstner & G. Riidiger: Local simulations of the turbulent ISM 



(200 5J). Earlier app roaches ([Brandenburg & Sokolofl2002t 
Kowal et al .1120051) were based on least square fit methods. 
The major drawback with these was that in regions where 
B or VB vanishes the inversion becomes singular. This can 
be circumvented by solving equation (Q~|) for fixed test fields 
with simple, well behaved geometry and gradients. To 
obtain the related EMFs one has to evolve an extra (passive) 
induction equation for the associated fluctuations: 

dtB[ v) - Vx [u'xB {v) + {u+qSlxy)xB[ v) 



u'xB[ u) - V VxB{ u) ] 



V-B? 



0. 



(2) 



We implemented these additional equations within NIR- 
VANA employing the constrained transport paradigm to ex- 
actly satisfy the solenoidal constraint. The actual method 
uses up-winding to guarantee stability while second order 
in space is attained via piecewise linear reconstruction. For 
this we apply the same slope limiter as in the actual code. 
Our procedure is very similar to t he methods described in 
Teyssier. Fromang & Dormyl ( 20061) . 

For the particular choice of the four t est fields Br v ) 
we use the ones from iBrandenburd ([2005), i.e., the low- 
est Fourier modes in the vertical direction. For each of the 
test fields we compute the corresponding mean electromo- 
tive force 

S {v) = u'xB 



(3) 



The dynamo coefficients can then be computed via equa- 
tion ((TJ to which the solution can be compactly written as 



&l»7ij3 



cos(fciz) sin(fciz 
-sin(fcjz) cos(fciz) 




(2j-2) 
(2J-1) 



, (4) 



with i,j G {1, 2}. In contrast to the least square fit method 
equation © can be directly computed for each z, yielding 
vertical profiles for the dynamo parameters. 




Fig. 2 Phase space distribution of the SN-heated plasma. 
Adjacent plots show the mass- (thick) and volume- (thin 
line) distribution of density (upper panel) and pressure 
(right panel). We also plot isothermal contours (labeled 
inK), and the equilibrium cooling curve (dashed). 



While single SNRs are largely confined to the midplane, 
super-bubbles break out of the central disk and drive mod- 
erate vertical flows. Their dense shells that are further com- 
pressed by shocks will form clouds that can efficiently cool 
and will in turn rain back into the gravitational potential 
thus f orming what is termed the galactic fountain (Bregman 
19801) . The Mach number of the flow is mildly supersonic 



below 1 kpc and becomes subsonic for the outer parts. 

If we turn off the SN-clustering the morphology changes 
quite drastically. Instead of well confined SBs we see more 
disrupted features and chimney-like structures that channel 
strong vertical outflows. The velocity dispersion in the hot 
phase is twice as high as in the clustered case. These differ- 
ences demonstrate the importance of a proper modeling of 
clustered explosions. 



4 Results and discussion 



4.1 Thermal distribution and velocity dispersions 



Since our vertical stratification is in radiative equilibrium, 
we do not observe an initial collapse of the disk. Turbulence 
builds up smoothly and after lOOMyr reaches a quasi steady 
state. Fig. Q] shows vertical and horizontal slices through 
the simulation box at t = 112 Myr. Most of the material 
is contained in cold clumps forming a 150 — 200 pc wide 
disk. Close to the midplane the network of clumps and fil- 
aments is permeated by strong shocks from the SNe which 
are continously creating hot cavities. Looking at the lower 
panels (a) and (e) of Fig. [T]one can see that for the region 
around the midplane there exists a significant correlation 
between the density and the magnetic field amplitude. The 
inferred slope of the correlation is somewhat steeper than 
the Chandrasekhar- Fermi value of 0.5 but roughly consis- 
tent with |B| ~ p 2 / 3 as indicative of compressional ampli- 
fication. 



In Fig. |2] we show the distribution of the SN-heated plasma 
as a function of density and thermal pressure at a time t 
1 1 2 Myr. The two stable branches of the equilibrium curve 
are densely populated but there also exists a considerable 
amount of gas in the radiatively unstable regimes. 

Averaged velocity dispersions for the ISM phases are 
5kms -1 (cold), llkms -1 (cool), 25kms _1 (warm), and 
40 — 60kms _1 (hot). While t he values for the latter are 
consis tent with the findings of de Avillez & Breitschwerdtl 
(120051) . for the cold phases we fall short by a factor of two. 
We consider this a resolution issue but it might also be re- 
lated to the different form of the driving. As a function of 
galactic height the time averaged rms-velocity rises steeply 
from its midplane value of 20kms~ 1 to its peak value of 
55kms _1 which is reached at about 1 kpc and then falls 
off to about 30 km s _1 at 2 kpc. Assuming that diamagnetic 
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pumping follows the inverse gradient in turbulence inten- 
sity this implies an inward transport of magnetic flux for the 
central part of the disk. 

4.2 Dynamo parameters 

In Fig. [3] we plot the main a- and 77-coefficients for our stan- 
dard run (mod. STD). To our knowledge this is the first 
time such profiles have been obtained from direct simula- 
tions of SN-turbulence. The corresponding vertically aver- 
aged amplitudes are listed in Table [2]- for comparison we 
also list the peak values from FER98 (where r\ t is distin- 
guished into horizontal and vertical part). Dynamo num- 
bers C a — aH/rj t and Cq — VlH 2 jr\ t are computed with 
H = 0.8 kpc and r\ t values from the inner part of the disk. 
In accordance with analytical estimations based on these 
numbers, the runs without shear are still sub-critical, i.e., no 
field-amplification is observed. The a 2 £l-dynamo, however, 
does exponentially grow at a time scale of ~ 250 Myr. 

During the course of the present simulations we do not 
yet reach the equipartition field strength. This implies that 
the values for the dynamo parameters (as well as the pitch 
angle) are representative of the unquenched regime. Our 
values are still affected by fluctuations and a longer time- 
base is needed for a more accurate determination. However, 
we find some interesting trends in the data at hand: The an- 
tisymmetric part -f z of the a-tensor is negative (positive) 
in the northern (southern) hemisphere, i.e., pumping is di- 
rected inward. This is contrary to the predictions by Ferriere 
for the case of non-interacting SNRs. In our simulations the 
inward pumping is opposed by an outward advection of the 
field with the mean flow. As can be seen from the middle 
panel of Fig. |3]both terms have about the same magnitude 
for z < 1 kpc while in the corona the wind is dominating. 
The balance in the inner region could be seen as an indica- 
tion that the effects of turbulent pumping and advection are 
naturally linked. In consequence vertical transport processes 
might be of lesser importance than formerly believed. 

The importance of modeling spatially coherent SBs can 
be seen from comparing the first two rows of Table[2] In the 
case without clustered SNe (mod. NCL) we observe strong 
vertical streaming motions reflected in high values for ^ z 
and u z . Also the turbulent diffusivity is high in the disk 
midplane which is not the case for the other runs where r\ t 
increases with galactic height. From the model with type- 
II SNe only (mod. SN2) one can learn that the turbulent 
diffusivity predominantly arises from the more broadly dis- 
tributed type-I SNe - despite their lower rate by a factor 
of eight. The overall level of turbulent diffusion is about 
3-6 x 10 26 cm 2 s _1 . The off-diagonal components of the 77- 
tensor are somewhat s mal l er and cannot yet be determined 
accurately. Following |jin| d 19961) we crudely estimate that 
the measured amount of diffusion suffices to damp short 
wavelength MRI-modes for reasonably high (3p - a defi- 
nite conclusion, however, requires further investigations by 
means of combined direct simulations. 
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Fig. 3 Dynamo a- and 77-coefficients for the STD model. 
Quantities indicated by the ordinate labels are plotted 
in dark (aaR,...) respectively light (a^, ...) colors. 
Shaded areas indicate the rms-fluctuation when applying the 
method to 4 temporal sub-intervals. 



By comparing 



with the SOCA results by 



Riidi ger & Kichatinovl (|1_993) we can estimate the coher- 
ence time t and the Coriolis number fl* = 2tQ of our 
turbulent flow. We find a value of r w 3.6 Myr which is 
about a factor of three smaller than commonly assumed. The 
even lower value for our model NCL indicates that a higher 
coherence time might be achieved by a more realistic pre- 
scription for the modeling of SBs. In general the obtained 
profiles are consistent with their SOCA counterparts. Our 
obtained value for 7 = \aj>p.\/\ari<b \ agrees with the SOCA 
value of 7 « 2.5 in iRiidiger & Hollerbach (2004) and is 



smaller than the value of 7 « 6 in FER98. 



In the lower panel of Fig. [3] we overplot the scalar ex- 
pression for the turbulent diffusivity which matches the pro- 
files obtained from the simulations. The peak of the veloc- 
ity dispersion close to the midplane is an artifact due to the 
static distribution of the SNe resulting in a disruption of the 
central disk. This issue has meanwhile been resolved. 
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Table 2 Mean values of extracted parameters averaged over « 150 Myr. Values for r\ t apply to the inner (\z\ < 0.8 kpc) 
respectively outer region of the disk. Coherence time r and Coriolis number ft* are estimated from a comparison with 
SOCA-profiles for a^. The coefficients could not be obtained for model SHR. 



4.3 Pitch angles 

In our models without shear we observe pitch angles rang- 
ing from +10° in the far outer regions to —45° near the 
midplane. This is hardly surprising since also the radial and 
azimuthal dynamo effects are found to be of equal strength 
in this case. If we include galactic shear with q = — 1, pitch 
angles become solely negative reaching values up to —30° 
as consistent with observations. The minimum of —10° is 
found at the midplane. In comparison, Hanasz et al. (2006) 
in their simulations of a cosmic ray driven galactic dynamo 
find the azimuthal field to be dominating by a factor of ten 
which corresponds to pitch angles of only ~ 5°. 

5 Summary 

We have performed resistive MHD simulations of the dif- 
ferentially rotating, stratified local interstellar medium. By 
integrating a radiatively stable initial stratification we were 
able to avoid the disk collapse occurring in previous isother- 
mal models. We have further demonstrated that SNe can in- 
deed produce an a-effect that is not dominated by vertical 
pumping. Our key findings are as follows: 

- The turbulent velocity dispersion in our box increases 
with galactic height, i.e., the diamagnetic pumping is di- 
rected inward. 

- We observe a mean vertical outflow of the same mag- 
nitude as the turbulent pumping. In consequence this 
drastically reduces the effective vertical transport of the 
mean magnetic field. 

- The ratio 7 of vertical pumping over the a-effect is di- 
minished if type-II SNe are modeled as clustered events. 

- The radial pitch angles found in our simulations are con- 
sistent with observations and suggest an a 2 £!-dynamo. 
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